#include "StVenKir.hpp"

// variational formulation language
#include <life/lifevf/vf.hpp>

namespace Life
{


  template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
  void
  StVenantKirchhoff<Dim, Order, Entity>::updateResidual( const vector_ptrtype& X, vector_ptrtype& R )
  {

    using namespace Life::vf;

    mesh_ptrtype mesh = M_Xh->mesh();

    element_struct_type u( M_Xh, "u" ); u = *X;
    element_struct_type v( M_Xh, "v" );
    element_struct_type force( M_Xh, "f" );


    //    AUTO( eta, 1000*Px()*( Px() -4 ) /**(Px()-2.5)*/  ) ;
    //AUTO( eta, 1000*(Px()-0.4)*( Px() -0.45 )*chi( Px()>=0.4 && Px()<=0.45 ) ) ;
    //AUTO( eta, 0.2*(Px()-0.4)*( Px() -0.45 )*chi( Px()>=0.4 && Px()<=0.45 ) ) ;
    //AUTO( eta, 0.2*(Px()-1.4)*( Px() -1.45 )*chi( Px()>=1.4 && Px()<=1.45 ) ) ;
    AUTO( eta, 0.3*sin(M_PI*Px()/4.) ) ;


    force = vf::project(M_Xh, elements(mesh), eta*oneY());
    //Matrice Identitee
    AUTO( Id , oneX()*trans(oneX()) + oneY()*trans(oneY()) );

    AUTO( Fv , Id + gradv(u) );
    AUTO( F , Id + grad(u) );
    AUTO( Fvf , Id + gradv(force) ); 

    AUTO( Evf, 0.5*(gradv(force)+trans(gradv(force)) ) + 0.5*trans(gradv(force))*gradv(force) ); 
    AUTO( Ev, 0.5*(gradv(u)+trans(gradv(u)) ) + 0.5*trans(gradv(u))*gradv(u) );
    AUTO( E, 0.5*(grad(u)+trans(grad(u)) ) + 0.5*trans(grad(u))*grad(u) );

    AUTO( Svf, lambda*trace(Evf)*Id + 2*mu*Evf ); 
    AUTO( Sv, lambda*trace(Ev)*Id + 2*mu*Ev );
    AUTO( S, lambda*trace(E)*Id + 2*mu*E );

    AUTO( dFf, gradv(force)); 
    AUTO( dF, gradv(u));
    AUTO( dFF, grad(u));

    AUTO( dEf, 0.5*(trans(Fvf)*dFf + trans(dFf)*Fvf) ); 
    AUTO( dE, 0.5*(trans(Fv)*dF + trans(dF)*Fv) );
    AUTO( dEE, 0.5*(trans(Fv)*dFF + trans(dFF)*Fv) );

    AUTO( dSf, lambda*trace(dEf)*Id + 2*mu*dEf ); 
    AUTO( dS, lambda*trace(dE)*Id + 2*mu*dE );
    AUTO( dSS, lambda*trace(dEE)*Id + 2*mu*dEE );
    /*
    //force volumique :
    AUTO(f,vec(cst(0.),cst(0.)) );
    */
    //--------------------------------------------------------------------------------------------------//

    form1( M_Xh, R, _init=true ) =
      integrate( elements(mesh), _Q<4*(Order)>(),
		 trace((Fv*Sv)*trans(grad(v))) );

    double penalisation_bc=this->application()->vm()["bccoeff"].template as<double>();
#if 0
    form1( M_Xh, R ) +=
        //integrate( elements(mesh), _Q<2*Order>(),
        integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<2*Order>(),
                   - trans(eta*oneY())*id(v) );
#endif

#if 1
    form1( M_Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Fixe")), _Q<4*Order>(),
                   //terme classique
                   - trans(Fv*Sv*N())*id(v)
                   //terme de symetrisation de la jacobienne
                   - trans((dFF*Sv + Fv*dSS)*N())*idv(u)
                   //- trans((dF*Sv + Fv*dS)*N())*id(v)
                   // - trans(F*S*N())*idv(u)
                   + penalisation_bc*trans(idv(u))*id(v)/hFace() );
#endif
#if 1
    form1( M_Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                   //terme classique
                   - trans(Fv*Sv*N())*id(v)
                   //terme de symetrisation
                   - trans((dFF*Sv + Fv*dSS)*N())*idv(u)
                   //- trans((dF*Sv + Fv*dS)*N())*id(v)
                   //- trans(F*S*N())*idv(u)
                   + penalisation_bc*trans(idv(u))*id(v)/hFace() );

    form1( M_Xh, R ) +=
        integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                   //terme de symetrisation
                   + trans((dFF*Sv + Fv*dSS)*N())*idv(force)
                   //+ trans((dFf*Svf + Fvf*dSf)*N())*id(v)
                   // - trans(F*S*N())*idv(u)
                   //!!!- penalisation_bc*trans(idv(force)/*idv(u)*/)*id(v)/hFace() );
                   - penalisation_bc*trans( idv(force)/*eta*oneY()*/)*id(v)/hFace() );
#endif

    double Newmark_gamma=0.5;
    double Newmark_beta=1/4.;
    //vf::pow(
    AUTO( cst_vel_displ , /*cst*/(Newmark_gamma/(Newmark_beta*M_bdf_displ_struct->timeStep())) );
    AUTO( cst_vel_vel ,  /*cst*/( (Newmark_gamma/Newmark_beta)-1 ) );
    AUTO( cst_vel_acc , /*cst*/( M_bdf_displ_struct->timeStep()*((Newmark_gamma/(2*Newmark_beta))-1 ) ) );

    element_struct_type buzz( M_Xh, "u" );
    buzz.add(cst_vel_displ,u);
    buzz.add(-cst_vel_displ, M_bdf_displ_struct->unknown(0));
    buzz.add(-cst_vel_vel, M_bdf_velocity_struct->unknown(0));
    buzz.add(-cst_vel_acc, M_bdf_acceleration_struct->unknown(0));

    AUTO( cst_acc_displ, 1./(Newmark_beta*std::pow(M_bdf_displ_struct->timeStep(),2)) );
    AUTO( cst_acc_vel, 1./(Newmark_beta*M_bdf_displ_struct->timeStep()) );
    AUTO( cst_acc_acc, ((1./(2*Newmark_beta))-1) );
    element_struct_type buzz2( M_Xh, "u" );

    buzz2.add(cst_acc_displ,u);
    buzz2.add(-cst_acc_displ, M_bdf_displ_struct->unknown(0));
    buzz2.add(-cst_acc_vel, M_bdf_velocity_struct->unknown(0));
    buzz2.add(-cst_acc_acc, M_bdf_acceleration_struct->unknown(0));

#if 0
    form1( M_Xh, R ) +=
      integrate( elements(mesh), _Q<2*Order>(),
                 trans(idv(buzz2))*id(v) );
#endif
//buzz = vf::project( M_Xh, elements(mesh), cst_vel_displ*(idv(u)-idv(M_bdf_displ_struct->unknown(0)) ) );

    //buzz->zeros();

    //form2( Xh, Xh, A ) +=
    // integrate( elements(mesh), _Q< 2*Order >(),
	//	 + trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0));


    R->close();

  }

} // Life


// instantiation
template class Life::StVenantKirchhoff<2,2,Life::Simplex>;
//template class Life::StVenantKirchhoff<3,2,Life::Simplex>;
